LXR irradiation scRNAseq Visium results
18 July, 2024
library(Seurat)
library(Matrix)
library(future)
library(future.apply)
library(ggplot2)
library(dplyr)
library(cowplot)
library(grid)
getGenePlots = function(gene, imgalpha = 0, qmax = NULL, topcolor = "#FF0000", lowcolor = "#FFFFFF", pt.alpha = 1){
g = SpatialFeaturePlot(DATA, gene, stroke = NA, combine = FALSE, image.alpha = imgalpha, alpha = pt.alpha)
if(!is.null(qmax)){
qval = quantile(do.call("rbind",lapply(g,function(x){x$data}))[,3],qmax)
g = SpatialFeaturePlot(DATA, gene, stroke = NA, combine = FALSE, image.alpha = imgalpha, max.cutoff = qval, alpha = pt.alpha)
}
if(imgalpha == 0){
g = lapply(g, function(x){x + theme(panel.background = element_rect(fill = "#CCCCCC"), panel.grid = element_blank())})
}
genemax = max(do.call("rbind",lapply(g,function(x){x$data}))[,3])
sfill = scale_fill_gradient(limits = c(0,genemax), low = lowcolor, high = topcolor)
return(lapply(g,function(x){x + sfill}))
}
plotGene = function(gene, imgalpha = 0, extratext = NULL,qmax = NULL, topcolor = "#FF0000", lowcolor = "#FFFFFF", pt.alpha = 1){
g = getGenePlots(gene = gene, imgalpha = imgalpha, qmax = qmax, topcolor = topcolor, lowcolor = lowcolor, pt.alpha = pt.alpha)
return(plot_grid(get_legend(g[[1]]), plot_grid(plotlist = c(lapply(g, function(x){x + NoLegend()})), ncol = 4), textGrob(extratext), ncol = 1, rel_heights = c(0.2,0.9,0.3)))
}
plotGeneColoc = function(gene1 , gene2, extratext = NULL,qmax = NULL){
g1 = SpatialFeaturePlot(DATA, gene1, stroke = NA, combine = FALSE, image.alpha = 0)
g1 = lapply(g1, function(x){colnames(x$data)[3] = gene1; return(x)})
if(!is.null(qmax)){
qval = quantile(do.call("rbind",lapply(g1,function(x){x$data}))[,gene1],qmax)
g1 = SpatialFeaturePlot(DATA, gene1, stroke = NA, combine = FALSE, image.alpha = 0, max.cutoff = qval)
g1 = lapply(g1, function(x){colnames(x$data)[3] = gene1; return(x)})
}
g2 = SpatialFeaturePlot(DATA, gene2, stroke = NA, combine = FALSE, image.alpha = 0)
g2 = lapply(g2, function(x){colnames(x$data)[3] = gene2; return(x)})
if(!is.null(qmax)){
qval = quantile(do.call("rbind",lapply(g2,function(x){x$data}))[,gene2],qmax)
g2 = SpatialFeaturePlot(DATA, gene2, stroke = NA, combine = FALSE, image.alpha = 0, max.cutoff = qval)
g2 = lapply(g2, function(x){colnames(x$data)[3] = gene2; return(x)})
}
colors = lapply(1:length(g1),function(x){rgb(g1[[x]]$data[,gene1]/max(g1[[x]]$data[,gene1]),g2[[x]]$data[,gene2]/max(g2[[x]]$data[,gene2]),0)})
g = lapply(1:length(g1), function(x){
#g1[[x]]$data[,3] = as.factor(rownames(g1[[x]]$data))
g1[[x]]$data$cell = as.factor(rownames(g1[[x]]$data))
names(colors[[x]]) = rownames(g1[[x]]$data)
g1[[x]]$mapping$fill = g1[[x]]$data$cell
g1[[x]] + scale_fill_manual(values = colors[[x]]) + NoLegend()
})
l1 = get_legend(ggplot(do.call("rbind",lapply(g1,function(x){x$data})),aes_string(fill = gene1)) +
geom_point(x = 1, y = 2) +
scale_fill_gradient(low = "black",high="red") + theme(legend.position = "bottom"))
l2 = get_legend(ggplot(do.call("rbind",lapply(g2,function(x){x$data})),aes_string(fill = gene2)) +
geom_point(x = 1, y = 2) +
scale_fill_gradient(low = "black",high="green") + theme(legend.position = "bottom"))
return(plot_grid(plot_grid(l1,l2),
plot_grid(plotlist = c(lapply(g, function(x){x + NoLegend()})), ncol = 4),
ncol = 1, rel_heights = c(0.2,0.9,0.3)))
}1 Load data
DATA = readRDS("Visium_dataset_irradiation.rds")
DATA$orig.ident = factor(DATA$orig.ident,
levels = c("STD_irrad0", "GW_irrad0",
"STD_irrad3", "GW_irrad3"))2 H&E images
To reduce file size, the H&E images have been removed from most plots. The images are shown here and can easily be combined with the spatial plots in e.g. Affinity.
g = SpatialFeaturePlot(DATA, "Abca1", stroke = NA, combine = FALSE, alpha = 0)
print(plot_grid(NULL,NULL, plot_grid(plotlist = c(lapply(g, function(x){x + NoLegend()})), ncol = 4), ggplot(x = 1, y = 1) , NULL, NULL, nrow = 3, ncol = 2, rel_widths = c(3,1), rel_heights = c(0.2,0.9,0.3)))
3 cNMF 3 factors
cNMF was performed using the Python package https://github.com/dylkot/cNMF. The results can be found at in the cNMF directory in the same directory as this report.
usages = read.table("cNMF/all_cNMF.usages.k_3.dt_0_5.consensus.txt", header = TRUE)
spectra = read.table("cNMF/all_cNMF.gene_spectra_score.k_3.dt_0_5.txt", header = TRUE)
usages_1 = apply(usages,1, function(x){x/sum(x)})
# create a new assay to store ADT information
#NMF_assay <- CreateAssayObject(counts = t(usages))
NMF_assay <- CreateAssayObject(counts = usages_1)
# add this assay to the previously created Seurat object
DATA[["NMF"]] <- NMF_assay
DATA@active.assay = "NMF"
factors <- rownames(DATA)
allfactors = list()
colors = c("#FF0000","#00FF00","#0000FF")
for(i in factors){
allfactors[[i]] = getGenePlots(i, topcolor = colors[as.numeric(gsub("X","",i))], qmax = 0.999, lowcolor = "#000000", imgalpha = 1, pt.alpha = 0.5)
}
plot_grid(plot_grid(plotlist =
lapply(1:length(allfactors[[1]]),function(j){
plotcolors = lapply(allfactors, function(x){
thisplot = x[[j]]
ggplot_build(thisplot)$data[[1]]["fill"]
})
colordf = do.call("cbind",plotcolors)
colors = apply(colordf, 1, function(crow){
#do.call("rgb",as.list(apply(col2rgb(crow)^2,1, mean)/255^2))
do.call("rgb",as.list(apply(col2rgb(crow)/255,1,max)))
})
g1 = allfactors[[1]][[j]]
#g1[[x]]$data[,3] = as.factor(rownames(g1[[x]]$data))
g1$data$cell = as.factor(rownames(g1$data))
names(colors) = rownames(g1$data)
g1$mapping$fill = g1$data$cell
g1 = g1 + scale_fill_manual(values = colors) + NoLegend()
return(g1)
}), nrow = 1
),
plot_grid(get_legend(allfactors[[1]][[1]]),get_legend(allfactors[[2]][[1]]),
get_legend(allfactors[[3]][[1]]), ncol = 3),
ncol = 1, rel_heights = c(3,1))
for(i in factors){
inum = gsub("X","",i)
topgenes = as.data.frame(t(sort(spectra[inum,], decreasing = TRUE)[1:15]))
colnames(topgenes) = "Loading"
topgenes$Gene = factor(rownames(topgenes), levels = rev(rownames(topgenes)))
tgbar = ggplot(topgenes, aes(x = Loading, y = Gene)) + geom_bar(stat = "identity") +
theme(panel.background = element_blank(), panel.grid = element_blank()) +
labs(title = paste("Factor",inum))
print(plot_grid(plotGene(i, qmax = 0.999),
plot_grid(tgbar,NULL, rel_heights = c(1.1,0.3), nrow = 2),
rel_widths = c(3,1)))
}
DATA = SetIdent(DATA, value = "orig.ident")
d0 = data.frame(p.value = rep(NA,3), diff = rep(NA,3))
d3 = data.frame(p.value = rep(NA,3), diff = rep(NA,3))
for(i in 1:3){
d0$p.value[i] = wilcox.test(DATA[["NMF"]]@counts[i,DATA$orig.ident=="GW_irrad0"],
DATA[["NMF"]]@counts[i,DATA$orig.ident=="STD_irrad0"])$p.value
d0$diff[i] = mean(DATA[["NMF"]]@counts[i,DATA$orig.ident=="GW_irrad0"]) -
mean(DATA[["NMF"]]@counts[i,DATA$orig.ident=="STD_irrad0"])
d3$p.value[i] = wilcox.test(DATA[["NMF"]]@counts[i,DATA$orig.ident=="GW_irrad3"],
DATA[["NMF"]]@counts[i,DATA$orig.ident=="STD_irrad3"])$p.value
d3$diff[i] = mean(DATA[["NMF"]]@counts[i,DATA$orig.ident=="GW_irrad3"]) -
mean(DATA[["NMF"]]@counts[i,DATA$orig.ident=="STD_irrad3"])
}
plot_grid(plotlist = c(lapply(VlnPlot(DATA,features = factors, pt.size = 0, combine = FALSE),
function(x){
x + NoLegend() +
geom_point(alpha = 0.1,
position = position_jitter(width = 0.3),
size = 0.1) +
geom_boxplot(width = 0.2, outlier.size = 0)}),
list(DotPlot(DATA, features = c("X1","X2","X3")) +
theme(legend.title = element_text(size = 8),
legend.text = element_text(size = 8))))
)
GW vs STD d0
| p-value | Difference in mean |
|---|---|
| 0e+00 | -0.0508736 |
| 0e+00 | 0.0231163 |
| 2e-07 | 0.0277573 |
GW vs STD d3
| p-value | Difference in mean |
|---|---|
| 0 | -0.0065889 |
| 0 | -0.0396660 |
| 0 | 0.0462549 |
topweights = as.data.frame(do.call("rbind",lapply(factors,function(i){
inum = gsub("X","",i)
df = as.data.frame(t(spectra[,names(sort(spectra[inum,], decreasing = TRUE)[1:20])]))
df$gene = rownames(df)
df
})))
topweights$gene = factor(topweights$gene, levels = topweights$gene)
heatmapdf = reshape2::melt(topweights)
colnames(heatmapdf) = c("Gene","Factor","Weight")
heatmapdf$Factor = factor(heatmapdf$Factor, levels = rev(c(1,2,3)))
ggplot(heatmapdf, aes(y = Factor, x = Gene)) +
geom_tile(aes(fill = Weight)) +
theme(axis.text = element_text(size = 8),
axis.title = element_text(size = 8),
legend.text = element_text(size = 8),
legend.title = element_text(size = 8),
axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, size = 6),
axis.ticks = element_blank(),
panel.grid = element_blank(),
legend.key.height = unit(10,"pt"),
legend.key.width = unit(5,"pt"),
panel.background = element_blank()) +
scale_fill_viridis_c(option = "magma")
sessionInfo()## R version 4.1.3 (2022-03-10)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: macOS Big Sur/Monterey 10.16
##
## Matrix products: default
## BLAS/LAPACK: /Users/jenfra/miniconda3/envs/das_nature_2024_2/lib/libopenblasp-r0.3.27.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] cowplot_1.1.1 dplyr_1.1.2 ggplot2_3.4.2
## [4] future.apply_1.11.0 future_1.32.0 Matrix_1.5-4.1
## [7] SeuratObject_4.1.3 Seurat_4.3.0
##
## loaded via a namespace (and not attached):
## [1] Rtsne_0.16 colorspace_2.1-0 deldir_1.0-9
## [4] ellipsis_0.3.2 ggridges_0.5.4 spatstat.data_3.0-1
## [7] farver_2.1.1 leiden_0.4.3 listenv_0.9.0
## [10] ggrepel_0.9.3 fansi_1.0.4 codetools_0.2-19
## [13] splines_4.1.3 cachem_1.0.8 knitr_1.43
## [16] polyclip_1.10-4 jsonlite_1.8.5 ica_1.0-3
## [19] cluster_2.1.4 png_0.1-8 uwot_0.1.14
## [22] shiny_1.7.4 sctransform_0.3.5 spatstat.sparse_3.0-1
## [25] compiler_4.1.3 httr_1.4.6 fastmap_1.1.1
## [28] lazyeval_0.2.2 cli_3.6.1 later_1.3.1
## [31] htmltools_0.5.5 tools_4.1.3 igraph_1.4.2
## [34] gtable_0.3.3 glue_1.6.2 RANN_2.6.1
## [37] reshape2_1.4.4 Rcpp_1.0.10 scattermore_1.1
## [40] jquerylib_0.1.4 vctrs_0.6.2 nlme_3.1-162
## [43] spatstat.explore_3.2-1 progressr_0.13.0 lmtest_0.9-40
## [46] spatstat.random_3.1-5 xfun_0.39 stringr_1.5.0
## [49] globals_0.16.2 mime_0.12 miniUI_0.1.1.1
## [52] lifecycle_1.0.3 irlba_2.3.5.1 goftest_1.2-3
## [55] MASS_7.3-58.3 zoo_1.8-12 scales_1.2.1
## [58] promises_1.2.0.1 spatstat.utils_3.0-3 parallel_4.1.3
## [61] RColorBrewer_1.1-3 yaml_2.3.7 reticulate_1.30
## [64] pbapply_1.7-0 gridExtra_2.3 sass_0.4.6
## [67] stringi_1.7.12 highr_0.10 rlang_1.1.1
## [70] pkgconfig_2.0.3 matrixStats_1.0.0 evaluate_0.21
## [73] lattice_0.21-8 ROCR_1.0-11 purrr_1.0.1
## [76] tensor_1.5 labeling_0.4.2 patchwork_1.1.2
## [79] htmlwidgets_1.6.2 tidyselect_1.2.0 parallelly_1.36.0
## [82] RcppAnnoy_0.0.20 plyr_1.8.8 magrittr_2.0.3
## [85] bookdown_0.34 R6_2.5.1 generics_0.1.3
## [88] withr_2.5.0 pillar_1.9.0 fitdistrplus_1.1-11
## [91] survival_3.5-5 abind_1.4-5 sp_1.6-1
## [94] tibble_3.2.1 KernSmooth_2.23-21 utf8_1.2.3
## [97] spatstat.geom_3.2-1 plotly_4.10.2 rmarkdown_2.22
## [100] data.table_1.14.8 digest_0.6.31 xtable_1.8-4
## [103] tidyr_1.3.0 httpuv_1.6.11 munsell_0.5.0
## [106] viridisLite_0.4.1 bslib_0.5.0